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We investigate the response properties of granular media in the framework of the so-caUed 
Random Tetris Model. We monitor, for different driving procedures, several quantities: the evolution 
of the density and of the density profiles, the ageing properties through the two-times correlation 
functions and the two-times mean-square distance between the potential energies, the response 
function defined in terms of the difference in the potential energies of two replica driven in two 
slightly different ways. We focus in particular on the role played by the spatial inhomogeneities 
(structures) spontaneously emerging during the compaction process, the history of the sample and 
the driving procedure. It turns out that none of these ingredients can be neglected for the correct 
interpretation of the experimental or numerical data. We discuss the problem of the optimization of 
the compaction process and we comment on the validity of our results for the description of granular 
materials in a thermodynamic framework. 

PACS numbers: 45.70.Cc (Static sandpiles; granular compaction), 05. 10. -a (Computational methods in statistical 
physics and nonlinear dynamics), 05.70.Ln (Nonequilibrium and irreversible thermodynamics) 

I. INTRODUCTION 

Granular media |l] are usually considered non-thermal systems because their thermal energy is so negligibly small 
with respect to other energy contributions (e.g. potential energy) that for all the practical purposes they are virtually 
at zero temperature. This feature draws a lot of consequences from the point of view of the validity of thermodynamics 
for such systems. One of the most important consequences is that, unless perturbed in some way (e.g. driving energy 
into the system), a granular system cannot explore spontaneously its phase space but it remains trapped in one of the 
numerous metastable configurations. One has then to look at the dynamics of a granular system always as a response 
to some perturbations and in general the response will depend in a non-trivial way on the rheological properties of the 
medium, on the boundaries, on the driving procedure and, last but not the least, on the past history of the system. 

In this paper we focus on the response properties of a class of lattice of models, the so-called Random Tetris Model 
(RTM) that, despite their apparent simplicity are able to reproduce many features of real granular materials: 
slow- relaxations during compaction |^,|| , segregation [Q , dilatancy properties , ageing Q . 

We shall be concerned in particular by the interplay between the response properties and the spatial structures 
that spontaneously emerge as a consequence of the dynamics imposed to these systems. Several examples in this 
direction have already been drawn: the phenomenon of structures formation that plays a crucial role in explaining 
the complex features of internal avalanching |Q or the coarsening phenomena that parallel the vibration-induced 
compaction process (to be discussed in a forthcoming paper 1^). 

We shall consider several procedures of vibro-compaction and we shall study how the system responds to different 
procedures and how the emerging inhomogeneities affect the response properties of these systems. The vibration 
procedures are specified in terms of the temporal function describing the evolution of the shaking amplitude. The 
simplest case is the one where one keeps the shaking amplitude indefinitely constant. More generally we shall consider 
complex procedures corresponding to sequences of cooling and annealing processes. In all the different cases we monitor 
several quantities. On the one hand, we focus our attention on global quantities as the global density, the response 
and the correlation functions. On the other hand we monitor some local quantities that allow us to investigate the 
large-scale structures spontaneously emerging in these systems as a response to the imposed perturbation (driving). 
The comparison between global and local quantities will be a valuable tool towards an understanding of how granular 
materials respond to perturbations and in this perspective of the importance of spatial structures. 

In this framework we can address several questions. One of the first issues we can investigate is the importance of 
the history, i.e. the specific procedure undergone by the sample before we perform our measurements. It is interesting 
to ask where the history of the system is encoded and, for instance, whether the global density represents a good 
parameter for the description of a static packing or one needs to specify other parameters. The analysis of the 
correlation functions will allow us to investigate under which conditions the systems exhibit ageing behaviour. On the 
other hand, with the knowledge of the effect of different kinds of perturbation, we shall address the problem of the 
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optimization of the compaction procedures that can be stated as follows: single out the best sequence of perturbations 
to impose to the system in order to maximize the density measured in a suitable part of the system. 

Our results will also allow us to comment on the existence of transitions in the response properties and on recently 
published results concerning the violation of the fluctuation-dissipation theorem (FDT) for granular media . 

The outline of the paper is as follows. In sect. II we recall the model definition, we define the different dynamical 
procedures and we introduce the quantities we shall look at for the analysis of the response properties. Sect. Ill and 
IV are devoted to the analysis of the response properties during a process of continuous shaking at constant amplitude 
and during cyclic procedures respectively. In sect. V we discuss globally all the results comparing different procedures 
and different histories and commenting on the consequences from the thermodynamic point of view. Finally in sect. 
VI we draw the conclusions. 



II. MODEL DEFINITION 

The essential ingredient of the RTM 1^ is the geometrical frustration that for instance in granular packings is due 
to excluded volume effects arising from different shapes of the particles. This geometrical feature is captured in this 
class of lattice models where all the basic properties are brought by the particles and no assumptions are made on the 
environment (lattice) . The interactions are not spatially quenched but are determined in a self-consistent way by the 
local arrangements of the particles. Despite the simplicity of their definition, these systems present an highly complex 
phase-space and their dynamics generates automatically a very rich gallery of time-space correlations: time-scales, 
spatial structures, memory, etc.. Furthermore they show a very interesting interplay between the dynamics and the 
time-space structures. It is worth to notice how in this class of models the origins of the randomness and of the 
frustration coincide because both are given in terms of the particle properties. 

Let us recall briefly the definition of the model, which includes, like in the real computer game Tetris, a rich variety 
of shapes and sizes. On a lattice each particle can be schematized in general as a cross with 4 arms (in general the 
number of arms is equals to the coordination number of the lattice) of different lengths, denoted by Ine, InWi ^SE, 
Isw, chosen in a random way. An example of particle configuration on a tilted square lattice is shown in Fig. n^. 




FIG. 1. Sketch of a local arrangement of particles in the Random Tetris Model: each particles can be schematized in general 
as a cross with 4 arms of different lengths, denoted by Ine, Inw, Ise, Isw, chosen in a random way. 

The static properties of the system are then completely characterized giving the random numbers defining the 
particles. These numbers are given once for all and it is clear how in this way one has a complete freedom in the 
choice of the system and the models used in |3|j^ are particular cases of the general model where one has chosen the 
particle sizes and shapes in a deterministic way. 

The interactions among the particles obey to the general rule that one cannot have superpositions. For instance 
one has to check that for two nearest neighbours particles the sum of the arms oriented along the bond connecting 
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the two particles is smaller than the bond length. It turns out that in this way the interactions between the particles 
are not fixed once for all but they depend on the complexity of the spatial configuration. We shall come back on this 
point later on in connection with the interplay between the dynamics and the emergence of spatial structures. 

The extreme generality of the model definition allows a large variety of choices for the particles. Just to give the 
idea of how the system can be chosen let us consider two different versions of the RTM that we shall denote A and B: 

(A) system with random elongated particles. This case represents the direct generalization of the system considered 
in Instead of considering two types of elongated particles with fixed sizes we consider elongated particles 
whose size is chosen randomly according to the expressions Ine = h'W, Inw = hs, Ine = f ^ ± O.lcZjy, 
Inw = d — Ine and equivalently for the other orientation, i] represents a random variable uniformly distributed 
in [0,1]. 

(B) system with random particles with spherical symmetry. In this case we particles whose size is chosen randomly 
according to the expressions I^e = IsWi Inw — ^se, Ine — 5d±0.2d?7, l^w = ±0.2(11]. 

In all our simulations we have always used, without loss of generality, systems with particles of type A. With the 
above given rules one can define the allowed configurations. For instance introducing gravity one can decide whether 
a certain configuration (packing) is mechanically stable or not. For instance Fig. ^ shows examples of mechanically 
stable configurations under gravity for the cases A and B. In these cases a tilted square lattice was used to implement 
the existence of a preferential direction in the system. 




FIG. 2. Example of a stable packing configuration for a system of type A (left) and B (right) (see text). In all our simulations 
we used particles of type A. 

The dynamics will then consist in general in a diffusion constrained by the particles geometry. In the following 
sections we shall be concerned with the case of the vibration-induced compaction phenomenon . 

One of the possible ways to implement a vibration or a shaking procedure in the RTM model is to consider a 
Monte-Carlo dynamics where the particles diffuse on the lattice according to some rules, as explained below. One can 
actually choose the specific dynamics in several ways. In the following we discuss some examples. 

The system is initialized by filling the container idealized as a lattice with lateral periodic boundary conditions and 
a closed boundary at the bottom. The procedure of filling consists in inserting the grains at the top of the system, 
one at a time, and let them fall down, performing, under the effect of gravity, an oriented random walk on the lattice, 
until they reach a stable position, say a position in which they cannot fall further. This filling procedure is realized 
by the addition of one particle at a time and stops when no particle can enter the box from the top anymore. 

In the case of shaking the dynamics can be divided in several alternating steps where the system is perturbed by 
allowing the grains to move in any allowed directions with a probability pup to move upwards (with < pup < 1-0) and 
a probability pdown = 1 — Pup to move downwards. Each step lasts until a fixed number of N moves per particle have 
been attempted with a fixed value of a; = Pup/Pdown- The quantity x can be related to the adimensional acceleration 
F used in compaction experiments [|o) through the relation F ~ l/log{\/ ^/x). 

More precisely the single dynamical step consists of the following operations: 1) extracting with uniform probability 
a grain; 2) extracting a possible movement for this grain among the nearest neighbours according to the probabilities 
Pup and Pdown] 3) move the grain if all the possible geometrical constraints with the neighbours are satisfied. 
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One is then free to choose the desired sequence of steps depending on the quantities one wants to monitor or on 
the experimental procedures one wants to reproduce. 

The easiest possibihty is the continuous shaking reahzed by letting the system evolve continuously with a fixed and 
constant x. We shall discuss this case in Sect. III. Another possible procedure, explored in mimics the process of 
tapping and it is composed by two alternating steps. First, in a "heating" process (tapping) the system is perturbed 
by allowing the grains to move in any allowed directions with a probability Pup to move upwards and a probability 
Pdown = 1 — Pup to move downwards. After each tapping has been completed (i.e. a fixed number of N moves per 
particle have been attempted with a fixed value of x = Pup / Pdown) we allow the system to relax setting p^p = 0. The 
relaxation process ("cooling") is supposed to be completed when no particles can move anymore under the effect of 
gravity, i.e. unless Pup is switched on. After this relaxation the system is in a stable static state and one restarts the 
cycle. 

More generally one can define complex cycles where one changes the value of x during the evolution of the system 
according to some specific temporal function. This is the most general case that allows to impose to the system 
specific histories. We shall discuss this case in Sect. IV. 

Before describing in a detailed way the results let us define the different quantities we shall be monitoring during the 
shaking procedure. Denoting with («,j) the coordinates of a generic site, where i indicates the horizontal coordinate 
and j the vertical one along the direction of gravity, we can define with m{i,j) the mass content of the site in 
such a way that m{i,j) = 1 if the site contains a particle and m{i,j) = otherwise. With these definitions in mind 
we shall monitor the evolution of the following quantities: 

Density Profile The density profile gives the value of the density (averaged over horizontal layers) as a function 
of the height. It represents the simplest, though rough, way to characterize the inhomogeneities in the system: 
since the gravity is acting in the vertical direction, we concentrate on the heterogeneities that can occur in this 
direction. In formulas we have 

1 ^ 

pU,t) ^ -J2^^At) forje[0:Af] (1) 

i=l 

where L is the width of the system and M its maximal height; 

Average density The density of the packing, i.e. the fraction of sites occupied with respect to the total number of 
sites, is measured after each relaxation step and, in correspondence with real experiments, we plot the behaviour 
of this density as a function of time. In order to avoid finite-size effects we considered systems with a linear size 
of at least L = 50 sites and, in order to be sure to observe bulk effects, we measured the density in the lower 
25% or 50% of the system; 

Response function The response function is defined as the change in the potential energy of the system for a 
small change in the value of x. What we do in practice is evolve the system for a certain time t^. At tw we 
define a replica of the system whose mass content will be described by m[ ^ and we let the original system evolve 
with the same x and the replica with x' = x + Sx where Sx is small enough to keep us in the linear response 
approximation. At times larger than tj^ we monitor the evolution of the potential energy for the two replicas 
defined as P{t + t^) = Yl,ij "^ij(^ + ^w) ■ (j + 1) and Pr{t + t^) = J^ij ^i.ji^ + ^w) ■ [j + !)• The response of 
the system is defined as 

where Npart is the total number of particles in the system. We also monitor the density profile of the replica, 
P^Ui^ + ^w), and the diff"erence Ap{j,t + t^) — p^{j,t + tyj) — p{j,t + tw) between the profiles. Note that 
can also be written as R{{t + t^^t.^j) — Ylj ^p{j, t + t^) ■ {j + 1) • L/ Npart- 

Correlation functions We have considered two kind of measures of the temporal correlations. The two times 
mass-mass correlation function, defined as: 



C(t + t^,t^) = — ^{mij{t + t^) ■mi^j{ty,)) 

-^^ part ■ ■ 



(3) 
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and the mean-square distance between the potential energies at times and t + tm, defined as: 

B{t + t^,t^,) = i '-^ ^ [h(t + U,) - h{t^)) (4) 

\ part part J 

where the overbar indicates the average over different realizations and h(t) indicates the height of the center of 
mass at time t. 

It is worth noticing that we decided to measure the two-times mass-mass correlation function via eq. ^ instead 
of the the two-times correlation function for the global density Q because, as it will be clear in the following, the 
definition of the global density is somehow arbitrary in these systems with a preferential direction imposed by gravity. 

We have used system sizes L x M of 60 x 60, 120 x 60 and 60 x 120 to ensure that finite-size effects were irrelevant. 
We let the system evolve for 10^ Monte-Carlo steps per particle after a preparation (waiting at constant x or other 
more complicated histories) during up to 10® Monte-Carlo steps per particle. Averaging was performed on up to 1000 
samples. 



III. RESPONSE PROPERTIES OF A SYSTEM SUBJECT TO CONTINUOUS SHAKING 

In this section we present the results of the simulations for a system subject to continuous shaking at constant x. 

The mean density, calculated in the lower 25% and 50% of the box, grows slowly, from an initial value po- This 
value is not universal but depends on the choice of the particles and on the fraction of the system where one performs 
the average. Fig.s || show the evolution of the density measured in the lower 25% and 50% of the system for different 
values of a;. It is interesting to notice that there exists an optimal value of x that allows to obtain the maximal density 
in the fraction of the system considered for the averaging procedure. At small times this optimal value depends on 
time due to the crossing of the different curves. This fact already suggests that, in order to compactify the system, 
the optimal strategy will not be to keep x constant, but rather to vary it. This type of behaviour was also recently 
observed in the parking lot model |l|]. However, in our case the mechanism responsible for these phenomena can be 
understood looking at the density profiles (see below), a concept which is absent in the parking model. Moreover if 
we compare the data for the bulk density computed as 25% or 50% of the bulk, we see that the curves differ and 
that the optimal value of x is around x ~ 0.6 in the first case, and around x ~ 0.8 in the other. This feature can be 
explained by noticing the combined effect of two factors: on the one hand the bulk compactifies better for larger x: for 
example, in the data of | [lO| , taken at the bottom of the sample, the density is an increasing function of the shaking 
amplitude; on the other hand, increasing x, the interface becomes broader and can affect the global bulk density 
effectively reducing it. As a result there will be an optimal x which for instance is smaller considering the measure of 
the density in the lower 50% of the system with respect to the measure in the lower 25% because one needs a larger x 
in order to extend the interface deeper and deeper. This phenomenology represents an indication that heterogeneities 
are quite relevant and that particular attention needs to be taken in defining the observed quantities. Our finding of 
an optimal shaking amplitude is in fact in agreement with previous experimental results . We recall for this that x 
and the adimcnsional acceleration T can be linked by the relation T ~ l/logil/ ^/x). In compaction experiments px| ] 
it was observed that the density, measured on the lowest part of the system, was an increasing function of the shaking 
amplitude. However, the authors were expecting this density to "decrease upon further acceleration increase" , since 
their experiments probed only the regime of relatively low shaking intensity. Moreover, such a decrease is expected 
to be more apparent in two-dimensional systems [p^ , which is the case of our model. 
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FIG. 3. Time-evolution of the bulk density (lower 25% (left) and 50% (right) of the system) for several shaking amplitudes. 
The insets show the density at time t = 10* as a function of x. Note in both cases the existence of an optimal value of x, that 
depends on the fraction of the system where one measures the average density. 

The growth law for the density has been shown by various authors to be, experimentally [ p^ and for many models 
1^,^, well described by the functional form: 

/,\ Poo — PO /^s 

P^'^ =P-- l + BHl + tlr) 

where poo, the asymptotic density, B and r are fit parameters depending on the shaking amplitude. Also in our case 
eq. ^ is satisfied and one can actually collapse all the curves in a unique logarithmic function. Fig. ^ shows the collapse 
of the density curves measured in the lower 50% of the system (one gets similar results using the curves of the density 
measured in the lower 25% of the system) for three different values oi x {x — 0.1, 0.6, 0.7) and for different values of 
tui (i.e. the measure of the density begins at instead of 0, which of course changes the parameters in eq. (||)): the 
functional form is tested by plotting the rescaled density {p{t) — po) / {B{poo — p{t))) versus t — t/r, and comparing it 
to ln(l + i). 




FIG. 4. Collapses of the rescaled density {p{t) — po) /{B{poo — p(i))) versus the rescaled time t — t/r, where B, poo and r 
are fit parameters, for a; = 0.1 (left), x — 0.6 (center) and x = 0.7 (right), and for densities measured after various t^. 

In all cases the functional form (j^) is quite well satisfied with r simply proportional to t^,, Poo being a bell-shaped 
function of x [ p^ and B exhibiting a complex dependence on t^, and x. Several remarks are in order. The fitting 
parameters depend on the fraction of the system where one measures the density. On the other hand the logarithmic 
fit starts to be valid after a transient time of the order of 50 — 500 iterations depending on x and on the fraction 
of the system used for the measurements. This behaviour is quite understandable. Before compaction can start in 
the bulk a weak decompaction is in order which is stronger for higher x and which extends at longer times if one is 
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measuring the density in a deeper region. One can actually notice the decompaction process in Fig. ^ especially for 
large values of x. Only after the decompaction process ends one can hope eq. |^ be valid. This is the reason why we 
have considered the collapse of the density curves obtained after a suitable waiting time . 

A different way of measuring the compaction is to look at the mean height, or potential energy P defined in section 
II. In this case, the whole system is considered. We see in Fig. ^ that, also in this case, there exists an optimal value 
of X for which one has the lowest position of the center of mass. Also in this case we find an agreement with the 
experimental findings of Ref . pO| where the density measured from the height of the system was an increasing function 
of the shaking amplitude, while the authors expected the presence of a peak "upon further acceleration increase" . 




FIG. 5. Potential energy (or mean height) of the system, for various x and measured after t„ — 10*. 

It is however already clear how the previous measurements, which are averages over the whole (or over an extensive 
part) of the system, cannot represent a comprehensive information about the system. For example, in the measure 
of the mean height the interface gives a large contribution while it is practically negligible if one is interested in bulk 
properties. Therefore before turning our attention to a more detailed description of the system let us now study the 
two-times correlation function C{t + t^,tw) which, as usual in ageing phenomena fl^ , gives a determination of the 
age of the system, since it depends on t and on i^,. 

In Fig.s H and we show the behaviour of C{t + t^,tw) for several values of x and . We observe for this model the 
typical ageing behaviour, with a first part for t <^ approaching a quasi-equilibrium curve where time-translation 
invariance is respected (i.e., for t ^ t^, C(t + tw^tw) approaches a curve depending only on t) a plateau, and, at 
i 3> iu,, a second decay, dependent on tw, corresponding to ageing. We note that this behaviour, very common in 
ageing phenomena ||l^, is more realistic than the one of the simplest Tetris model [H, where the plateau is at C = 1, 
so that only the second decay is observed p. As grows, the plateau is better defined and the second decay appears 
at longer times. 

The global properties of the correlation function depend smoothly on the shaking amplitude: as x grows, the 
correlation decays faster, as shown in Fig. However, the curves do not differ qualitatively and have similar shapes. 
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FIG. 6. Correlation function C{t + tw,tw) as a function of t, for x = 0.1 (left), x — 0.6 (center) and x = 0.7 (right), and 
various values of tw The ageing behaviour with a two-steps relaxation is evident. 
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FIG. 7. Correlation function C{t + tw,tw) as a function of t, for x — 0.2 and x = 0.4, tw — 10 and t^, = 10 : the correlation 
decays faster for higher x. 

In j|] the two-times correlation function for the global density was measured in the simplest version of the Tetris 
model. In this case it was proposed that the relaxation was of the form ln(t)/ ln(tu,) (the exact proposed form was 
(1 — Coo) ln((ito -I- ts)/T)/\Ti{{t + tw + ts)/T) + Coo, where Coo, tg and r were three fit parameters), for very small values 
of X (in the range [10^^; 10^^]). This form could fit all the curves since, as already mentioned, the first decay was 
inexistent. 

Since, as it has been pointed out in the above discussion, the definition of the system density is somehow arbitrary, 
we decided to measure the two-times mass-mass correlation function via eq. p|. In our case, we only attempt to fit the 
second decay, i.e. the ageing part; since we expect weak ergodicity breaking ]16|, we propose a form going to zero at 
long times: 

C(i + i.,U.^-^^^-^-^fori»t. (6) 

with a, b, T fit parameters. We show in Fig.s |^ and |^ that we can collapse the curves using eq. ^ and plotting the 
rescaled correlation function C versus the rescaled time t according to 



C 



bC 



a+{Ab~l)C (^Y) 

i = t/r, 

and obtaining C = \/{A + ln(l + 1)). A is an irrelevant parameter introduced to avoid the divergence of (7 at f ^ 0, 
and whose value can be read in the figures. The behaviour of the fit parameters is the following: a is roughly constant, 
while r evolves proportionally to t^, and b as the inverse of the logarithm of This leads to the conclusion that 
the overall behaviour of C{t + t^^t^) is of the form log[tw + t)/log{ty^). We note that this is in agreement with the 
findings of [||, on a much wider range for x, but in contrast with the parking- lot model, for which a t/t^^ behaviour 
has been observed ||l^, where, as in Q, a different definition of the correlation function was used: the correlation 
function of the densities. An experimental measure of the correlation functions would be welcome to discriminate 
between these predictions. 
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FIG. 8. Collapse of the ageing part of the curves for the same values of x and (5 values for each x) as in Fig. |6| the 
rescaled correlation (see text) is plotted versus the rescaled time t/r, and compared to the form l/(^ + ln(l + 1)); the first part 
of the curves, not collapsed, corresponds to the first relaxation (approach to a quasi-equilibrium behaviour). 
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FIG. 9. As in Fig. g for = 0.2 and x = 0.4. 



For completeness, we show in Fig. 10 the behaviour of the mean-square distance between the potential energies, 
B{t + tyj, tyj) defined in the introduction though eq. ^. B{t + i^, tyj) is an increasing function of t, displaying ageing 
behaviour with two steps separated by a plateau, as for the correlation function. 
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FIG. 10. Evolution of the mean-square distance between the potential energies, B{t + tw,tjn), with t, for x — 0.2 (left), 
X — 0.4 (center) and x = 0.7 (right), and various tw 

Before turning our attention to the response function defined in eq. |2| let us add one important information con- 
cerning the density profiles. It turns out in fact that it is essential to look at the inhomogeneities in the system in 
order to correctly interpret the response results. We shall do this looking at the density profiles and at the differences 
between the profiles of the system and its perturbed replica, which will give information on the spatial structures, at 
least along the vertical direction. 

We have monitored the density profiles p{j, t), as defined in eq. |^, of the system and of its replica after i^j. While 
the curves of the evolution of the bulk density (or of other global quantities like C or B) have similar shapes for all 
values of x (Eq. ^ and Fig. the density profiles can exhibit very different behaviours. In Fig. ^l], we display the 
short-time {t — to t — 10"') evolution of the profile for several values of the shaking amplitude, starting from the 
same profile for all values of x. In Fig. |l^, we display the successive evolution, for < t < + 10^ with = 10^. 
The first observation is that the width of the interface is larger for higher shaking amplitudes, which is a quite intuitive 
result. Moreover, we see that, for small shaking amplitudes, a very dense layer forms just under the interface. This 
layer (whose height is of 5 to 10 lattice sites) is able to block the compaction process: in order to make the system 
more compact particles have to rearrange also in the bulk, and the dense layer acts as a barrier. The bottom part of 
the sample is therefore almost not evolving, as shown in Fig. [l^. As time evolves the dense layer becomes broader, 
though in a very slow way. In all cases, the comparison of Fig. |ll| and |2| also shows that the short-time dynamics is 
much faster than the successive evolution. 

We note that these results are in agreement with experimental results showing that, at shaking amplitude not too 
large, the compaction is more efficient in the higher parts of the media (see first reference in [|l^) and that the locally 
measured density is larger in higher parts of the sample. 
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FIG. 11. Density profiles for various values of x, evolving in time (t = to t = 10"'), with no waiting time: the evolution 
starts from the same profile for all values of x. We see the formation of the dense layer for small x, and the evolution towards 
a smoother profile for large x. 
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FIG. 12. Density profiles for various values of x, evolving in time (t = to t = 10^) after a waiting time = 10*. The 
interface is broader for higher x, and the dense layer below the interface is more pronounced and evolving more slowly for lower 
X. Comparing with Fig. Illl, we see that the overall evolution is much slower in all cases. 

For higher x, diffusion is easier, and the layer forming at the top is less dense, broader, and evolves faster (see in 
Fig. ^ the evolution of the dense layer for x = 0.1, 0.2, 0.4, 0.6); the compaction process in the bulk is therefore 
facilitated: thus, even if the interface is broader, which could indicate a less dense system, the bulk will be able to 
compactify much better. The global density can therefore be in fact enhanced. This explains that, at least at finite 
times, the curves of the density can be higher for higher x, while one would expect the contrary starting from a dense 
system; of course, for very strong shaking, the system is very loose and the interface is very wide: this explains the 
existence of an optimal shaking amplitude for compaction (see Fig. |^) and, as mentioned before, why this optimal 
amplitude depends on the way of defining the bulk (it roughly corresponds to the value of x for which the interface 
attains the defined bulk). 

We have again an indication of the importance of looking at the density profile: two very different profiles can have 
the same bulk density or potential energy (mean height). 

Let us now turn to the analysis of the response function. As described in the introduction, we let the system evolve 
at constant x during a certain waiting time t^,, and then make a copy of the system, which we submit to a slightly 
different shaking, i.e. x + dx. We have mostly used dx = 0.01, and checked with dx = 0.005 and dx = 0.02 that the 
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system was in the linear response regime (especially for x = 0.1 for which dx = 0.01 is 10% of x). 

The response function has been introduced in Sect. II through eq. ^: it is defined as the difference between the 
potential energy of the two copies of the system. Therefore, a positive response means that the perturbed system 
(at higher x) has a higher energy, and therefore it is less compact than the unperturbed one. This is what one 
clearly expect at high x for example, where the system is very loose and a higher x means that the particles are less 
subject to gravity. On the other hand, a negative response means that the perturbed system is more compact than 
the unperturbed one. This is not to exclude a priori, since we have already seen in Fig. ^ that the density at fixed 
time is not a monotonous function of x. 

For a; = 0.1, we see indeed (Fig. ^3|) that the response exhibits first a positive branch followed, after a certain time 
that depends on x and t^, by a negative branch. The system is first relatively decompactified by the perturbation 
(positive response), but then, at later times, it compactifies better. This phenomenon was first noted in the case of 
the Tetris model in Q , where however only small values of tm were used, and more emphasis was given to the negative 
part. For a fixed x, however, as t^ grows, the positive part extends to longer times and cannot be neglected. 




FIG. 13. Response function R{t + tw,tw) versus t, for several values of x and tw The vertical scale is the same for the left 
and right parts. The response tends to remains positive for a longer time higher are either x or 
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FIG. 14. Response function R{t + 1^,1^) versus t, for = 10'* (left) and = 10^ (right), showing that, at fixed t^, the 
response is an increasing function of x. 

As X grows, the positive part gets a higher amplitude and extends to longer times, and for x > 0.7 no negative 
response can be reached (Fig. ^ and |l|). 

This apparently odd behaviour can be understood by looking at the density profiles, and especially at the differences 
of the profiles between the perturbed and unperturbed systems, i.e. at the spatial distribution of the response. The 
response is made by two main contributions. On the one hand the interface gives a positive response, i.e. a larger x 
will naturally lead to a looser interface, because a larger x means higher probability for the particle to move upwards. 
On the other hand, however, the effect on the bulk is less obvious, since the particles can be blocked by the ones 
situated above (as is the case at small x, with the dense layer appearing), and need a global rearrangement of other 
particles in order to be able to move. 

We show in Fig. |l^ how the difference Ap(j, t + t^) — p^{j, t + t„,) — p(j, t + t^) as a function of j evolves in time 
after tu,. A positive Ap(j, i + t^,) means that the perturbation has locally (at height j) compactified better (giving 
a local positive contribution to the response function), while a negative Ap(j, t + t-^j) corresponds to a locally less 
compact perturbed system (local negative contribution to the response function). Besides, recalling the definition of 
the response function, R{t + t^,tw) — + tw)ij + ^)L/ Npart, it is important to remark that the values of 

Ap(j, t + tw) coming from the interface (large j) bring a strong contribution to the response due to the term [j + 1). 

At low x {x — 0.1), the bulk is blocked by the thick dense layer just below the interface, as shown in Fig. |l^. 
Therefore, the first effect of the perturbation dx is just to decompactify the interface, which gives rise to a positive R. 
However, once the interface is loosened, particles in the bulk may be allowed to rearrange more freely, and we obtain 
therefore a better compaction of the bulk, and a negative R. This phenomenon is shown in Fig. |l^ by the creation of 
a dip at heights j just below the interface, with a bump p^{j,t) > p{j,t) just above the interface at low t, and then 
by the decrease of this bump with a mass transfer toward the low-j values of the dip. Let us remark, however, that a 
negative response, which means a better compaction for the perturbed replica, does not give any substantial change 
in the profile: the response to the perturbation is very heterogeneous, and the bottom part of the sample simply does 
not feel it. As grows, the interface and the dense layer become more and more compact, and therefore harder and 
harder to decompactify. This explains why the response function stays positive for longer times: the decompaction 
process takes longer and longer (Fig. |l6| ). 
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FIG. 15. Temporal evolution {t — Otot — 10^ ) of the difference between the density profiles of the perturbed and unperturbed 
systems, Ap(j,t + t-uj), after t^j = 10* and various values of x. Comparing this figure with Fig. |l^ one observes that the dips 
and bumps are located near the interface, and they are larger for larger x, i.e. for a broader interface. 



Ap(j,t) 
0.08 
0.06 
0.04 
0.02 



x=0.1 



t = 10 



x=0.1 



t,„= 10 




20 40 



60 



20 40 60 80 



16 



FIG. 16. Temporal evolution (t = to t = tw + W ) of the difference between the density profiles of the perturbed and 
unperturbed systems, Ap(j, t + tm), for a; = 0.1 and two waiting times: tw = 10* and tw = 10®. The larger is tw, the slower is 
the overall processes. 

At higher values of x, the interface becomes smoother, the blocking layer becomes less dense and the effect of 
the perturbation on the interface becomes stronger and stronger. At strong shaking, i.e. large values of x, only the 
interface contributes to the response function, which therefore is positive at all times. We see in Fig. ^that the effect 
of the perturbation is to transfer particles upwards. 

In summary, as either tw or x grows, the bulk is more compact, and therefore its contribution to the response is 
smaller. The response function tends to become positive due to the contribution of the interface. 



IV. RESPONSE PROPERTIES FOR A CYCLIC SHAKING PROCEDURE 



In this section we consider a shaking procedure defined as a sequence of steps where x is varied following a cycle 
or a more general function. A generic procedure is defined giving an initial and a final value of x, xj and xpj the 
maximal value of x, Xmax, the value of the increments in x (in our simulations we have always used Ax — 0.01) and 
the time interval At the system spends at each value of x. We have considered values of At equal to 10^, 10'^, lO'^ 
and 10^. The two main procedures we have considered are defined in the following way: 

• cycle: increase from xi — 0.01 to x„iax = 0.8, and then decrease to xf- 

• cooling: simple decrease from xj = x^ax [xi — 0.6, 0.7, 0.8, 0.9) to xp. 

The final values of x were xp = 0.4, 0.2, 0.1. Once the final value xp of x was reached, x is kept constant and the 
measurements of the response function (with a copy evolving at a: + dx), the correlation function and the density 
profiles were done either without any further delay, or with an additional waiting time at constant x = xp of 10* or 
10^ time-steps. 

We shall see that such procedures, similar to a slow cooling for thermal systems, allow to reach large densities, 
seemingly unreachable at constant x. We note that these procedures are indeed the ones used experimentally to 
efficiently compactify granular systems. The examination of the density profiles will allow to gain a deeper insight for 
the effectiveness of such cycles. 

We first plot in Fig. ^ the evolution of the density during a cycle; the influence of At is clear, and similar to the 
case of cooling: larger At allow to reach higher densities. Moreover, the obtained densities are impressively higher 
than the densities obtained at constant shaking amplitude. 




FIG. 17. Bulk density versus x during the cycle (left) or cooling (right) procedures, for various values of Ar and Xma^- On 
the left figure is shown for comparison the behaviour of two cooling procedures with At = 10* and 10^. We see that a simple 
cooling from Xmai is equivalent to a cycle with the same Xma-j^ ■ In general higher densities are obtained with slower procedures 
(larger At) or higher values of 



17 



We also checked that the part of the cycle with growing x has no practical use or influence on the "cooling" part 
(see Fig. p^ . Only the maximal value of x is relevant. This seems reasonable, since structures formed at low x are 
destroyed by a shaking at larger x. On the other hand, if Ar is not very large, we have observed that a second cycle 
can be useful to obtain still higher densities. The right part of Fig. |l^ also shows that higher densities are obtained 
with higher values of Xmax ■ 

It is quite interesting to look at the density profiles at various xp after a cycle (see Fig. ^8|). For xp — 0.1, the bulk 
is much denser than after any reachable tw at constant x, and the interface is as steep. (Note however that, for At 
too small, there is still a dense layer at the interface for xp — 0.1.) Moreover, at fixed x/, a larger At yields better 
compaction in the bulk, as does, at fixed At, a larger x/. 

The comparison of the profiles at xp = 0.4, 0.2 and 0.1 (at fixed xi and At) shows that the bulk parts of the 
profiles are identical: only the interfaces change and they are steeper for lower x. When x is lowered, the bulk retains 
its properties while the interface is gradually sharpened. This means that, in order to better compactify, one has to 
take into account that high values of x are effective for the bulk while low values of x make the interface denser and 
steeper. Besides, the larger xj (see bottom pictures of Fig. |l^), the deeper the bulk is affected at the beginning of 
the cooling, and the more compact the system is at the end of the cooling. 



0.8 




FIG. 18. Top left: influence of xp on the density profile after a cycle, for Xmax = 0.8, Ar = 10*; the bulk is not affected, 
but the interface is steeper for smaller xf- Top right: at constant Xmax = 0.8 and xf ~ 0.1, effect of Ar; the bulk is denser 
for larger Ar, while the steepness of the interface is not changed; xf = 0.1 corresponds to a complete cycle, while x'j^ = 0.1 
corresponds to a cooling from Xmax ~ 0.8. We see that the two procedures are equivalent. Bottom left: influence of xi — Xmax, 
at constant Ar — 10*; the bulk is denser for higher a;/, the steepness of the interface is not changed. Bottom right: same for 
Ar = 10^ 
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From these observations we deduce that the optimal method to compactify is to begin with a high xi , and decrease 
X as slowly as possible. If Xi is higher, then one is allowed to take a lower value of At, thus gaining time in the 
compaction process. Note moreover that, for xj — 0.9, the difference between the profiles obtained at At = 10^ or 
10^ is very small. Moreover, if only the bulk has to be compactified, one can stop the process at intermediate values 
of X (see Fig. |l8|) while, in order to also have a sharp interface, the process has to be continued down to low values 
of X. It would be interesting to have detailed experimental tests of these predictions. 

After a cycle or a cooling procedure, the mass-mass two-times correlation function (Fig. 19) consistently shows an 
effective age of the system much larger than the real one: for example, for a global real time of 7 • 10^ (cooling from 
xi = 0.8 to xi = 0.1 in steps of Ax — 0.01 with At = 10"*), the effective age is much larger than 10^ (largest waiting 
time simulated at constant x). This fact is confirmed if we wait an additional = 10^ or 10^ after the cycle: we 
observe time-translation invariance, showing that the system is in an equilibrium or quasi-equilibrium state. 




FIG. 19. Correlation function after a cycle with xf = 0.1, Xmax ~ 0.8. The longer the cycle, the older the system looks 
(i.e. the slower the correlation decays). The curves for At = 10" with U = 0, 10*, 10^ are perfectly collapsed, showing 
time-translation invariance. 

The response function, shown in Fig. pOl, is also independent of t^, (for — 0, 10^, 10^), and is always positive. 
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FIG. 20. Response function after a cycle with xf ~ 0.1 and xf = 0.8, Xmax ~ 0.8, and a time = 0, lO**, 10^ at constant 
X = xf ~ 0.1. The response is small with large fluctuations, but definitely positive. We show only the case xf = 0.1, the 
response being an increasing function of a::F. 

Once again, this behaviour can be understood by the analysis of the difference of the density profiles of the two 
replicas. Fig. |2l| shows Ap{j,t) for several combinations of xj and xp. This figure is to be compared with Fig. |l^. 
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In all the cases one observes that the change in xp only modifies the interface, broadening it. The interface in this 
case is always composed by a denser layer (bump) lying above a looser layer (dip). The only effect of the perturbation 
is thus a broadening of the interface, even at small xp, thus causing a positive response, because the bulk has been 
already efficiently compactified by the cycling procedure. 
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FIG. 21. Temporal evolution (from t = to t = 10*') of the difference between the density profiles of the perturbed and 
unperturbed systems, Ap(j,t), after a cooling from xi — 0.8 to xp ~ 0.1, 0.2, 0.4, with Ar — lO'*, (top left, bottom left and 
bottom right respectively) or from xi = 0.9 to xp = 0.1 (with Ar = 10*) (top right). 



V. COMPARING DIFFERENT PROCEDURES: THE IMPORTANCE OF THE HISTORY 

In this section we try to give a unitary view for the results we have shown so far. Since in a granular system 
there is no dynamics without any energy injection, i.e. without any external perturbations, it is highly important 
to deeply understand the response properties to such perturbations. One of the most important point to stress 
is the fact that the response is never homogeneous. One can never assume that the properties of the system are 
homogeneously distributed, not even assuming a coarse grained point of view. The perturbations drive the system into 
an instability mechanism that generates large-scale spatial structures . Let us discuss what are the consequences 
of this phenomenon that has been named Self- Organized Structuring M. 

First of all, one is not allowed to describe a static packing in terms only of a scalar quantity like the density. It is 
evident how it is possible to construct several different packings corresponding to the same global average density, each 
one with completely different rheological properties (for a concrete example concerning the behaviour of a granular 
system subject to shearing see [^). Recent experiments [|l8| have made this point clear. Depending on the properties 
investigated one is then forced to enlarge the parameter space in order to give a reasonable description of the system. 
One important ingredient to take into account is the history, e.g. the ensemble of dynamical procedures the system 
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has undergone till the moment where we analyze it. One can then ask where the information about the past history 
is encoded and whether is possible to take into account this history in some suitable coarsed grained view. One 
possibility is to consider the density as a local parameter and describe the system in terms of some density map. It is 
precisely in this spirit that in this paper we have considered, as a first step in this direction, the analysis of the time 
evolution of the density profiles . This corresponds to investigating the properties of the system in the direction 
of the imposed external field: the gravity. 

The presence of large-scale structures in the system makes the problem of the response properties far from being 
trivial. One of the first question one could ask concerns the best procedure to compactify a given sample. The best 
strategy refers to the strategy that allows to obtain the highest density in a given time. Alternatively one could 
consider the best strategy to reach a given density in the smallest time. Since different parts of the sample (in the 
zero-th order schematization one can consider the bulk and the surface) respond in a different way to different values 
of the shaking amplitude x, it is obvious to conclude that the best strategy will not coincide with a procedure in 
which one keeps x fixed indefinitely. In this case, in fact, as we have shown in Sect. Ill, the average bulk density at 
fixed time is not a monotonous function of x. It is then quite intuitive to figure out different procedures where a; is a 
complex function of time. We have shown in Sect. IV that one of the best procedures is to consider a cooling process 
where one starts with a relatively high value of cc = Xmax and decrease progressively x with a rate 1/Ar up to a final 
value xf- The higher are either Xmax or At, the larger will be the final asymptotic density. The smaller is the value 
oi xf, the larger will be the compactified region and the sharper the interface. 

The rationale behind this definition of the optimal compaction procedure can be understood as follows. The best 
way to compactify globally the system is to start from the compaction of the bottom part of the system. In order 
to do this one should choose high values of x which allow to extend the interface of the system, i.e. the mobilized 
region, deeper and deeper. The procedure then proceeds reducing progressively the value of x in order to compactify 
sequentially regions at larger heights. It is clear that a small value of x cannot affect the part of the system already 
compactified with a larger value of x. That is why in order to better compactify also the system interface one has to 
continue the procedure to very small values of x. In this way one can associate to each particular value of x during the 
procedure the optimal compaction of a specific region of the system. On the other hand it is clear why the compaction 
process gives better results for larger At. Spending longer on a certain value of x allows to better compactify the 
region corresponding to this value of x. The shaking procedure at constant x cannot thus be effective because a large 
X will only be able to compactify the deep bulk while a small x will create a high density layer below the interface 
(see Fig. |l^) which will stop for a very long time the bulk compaction. 

The presence of large-scale structures draws a lot of consequences also in the behaviour of the response function as 
defined in eq. ^ From the results presented in the previous sections, it is evident as, even in a very simplified picture, 
the sign of the response function depends on a complex convolution of several contributions: the spatial structures 
spontaneously emerging in the system (again at the zero-th order of approximation one can consider the sum of two 
coupled contributions coming from the surface and the bulk), the value of the shaking amplitude and the past history 
of the system (encoded in the value of for a constant shaking amplitude or, more generally, in the whole definition 
of the dynamical procedure). It is then evident how trying to explain the results about the sign of the response 
function on the only basis of the shaking amplitude (x) can be fallacious. There does not exist, as suggested in 
a transition in x such that the the response function is positive above a certain value of x and negative below this 
value. In the case of constant shaking amplitude, no matter what the value of x is, one typically observes a positive 
response due to the contribution of the system surface, eventually followed by a negative regime at times that depend 
on X and tw- the larger are x or ty^,, the larger is the time over which one sees a positive response. This is due to 
the complex balance between the contributions coming from the surface and the ones from the bulk. These features 
were not considered in |^ where the role of inhomogeneities was neglected and where only small values of t^^ were 
considered. On the other hand if one looks at the response after a cyclic procedure, as described in sect. IV, one 
realizes that the response is always positive because in this case only the interface is giving an important contribution. 

It is interesting to compare the results obtained for the mass-mass two-times correlation function using different 
driving procedures. We have shown in Sect. Ill the presence of ageing in a system driven with a given constant 
X. In this case one observes a two-steps relaxation of the correlation function, typically observed in glassy systems. 
The second relaxation is quite well described by a function decaying as log(tw)/log(t + 1^]) no matter the value of x. 
Similar behaviour has been observed for the mean-square distance between the potential energies, B{t + ty^,tw)- On 
the other hand, the behaviour of the correlation function turns out to be completely different if the system is subject 
to a more complicate driving procedure. In Sect. IV we have shown that the use of a cyclic or of cooling procedure 
brings the system in a state where time-translation invariance holds, i.e. there is not anymore ageing. Experiments 
analyzing the presence or absence of ageing in granular material would certainly be welcome, and a first step is being 
taken in Q . In this case the state reached by the system is such that the bulk is almost completely decoupled from 
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the system interface, i.e. the bulk density does not change anymore with time. At this stage the dynamics of the 
system is concentrated on the interface which is almost in an equilibrium state compatible with the final value of x, 
xp. The system does not age anymore and the global response (see Fig. |2^) is always positive. For another approach 
to the relaxation dynamics of granular media that focuses on the asymptotic stationary state reached after a very 



long shaking procedure we refer the reader to |20 



These last observations arise a further question concerning the possibility to associate to a granular system a 
unique scalar temperature when the space translation invariance is broken. In pO[ | it has been shown that this is 
indeed possible if the system is in a stationary state. In this case one can show that it there exists a "temperature" -like 
quantity which is made uniform everywhere in the system and which is related do the derivative of a suitable free- 
energy-like functional. It is evident how one gets immediately into troubles if the system is far from being homogeneous 
and if the presence of spatial structures is accompanied by a breakdown of the time-translation invariance (ageing 
behaviour). The above discussion about the sign of the response function represents an example of how data analysis 
can be misleading. It is then clear how in general one cannot describe the response properties of a granular system 
in terms of a unique parameter, e.g. the temperature. Considering different values of x or or considering different 
dynamical procedures, one is simply exploring different non-equilibrium regions of the phase-space that correspond to 
different spatial inhomogeneities. From the above discussion it should be clear that a negative value of the response 
function does not mean that one can associate a negative effective temperature to the whole system (in the spirit of 
the Fluctuation-Dissipation Theorem (FDT) j2^ ) . In this case one has simply that for the specific values of x and 
explored, the contribution coming from the bulk is important. This contribution can change depending on the global 
history (which is encoded in the spatial structures), on x and and thus, unless one is able to introduce more subtle, 
and local, indicators, it seems to us quite meaningless and misleading to try to extend to this case the definition of 
an effective temperature. From this point of view, unless one is able to define suitable free-energy-like functionals p^ ] 
taking into account the heterogeneities of the system, any statement about the validity of the Fluctuation-Dissipation 
Theorem (FDT) in such an heterogeneous systems [^, statement which anyway has to be considered in the limit 
oo, seems to us hazardous. 

VI. CONCLUSIONS 

In this paper we have investigated the response properties of granular media in the framework of a recently proposed 
class of models, the so-called Random Tetris Model. On the one hand, we focused our attention on global quantities 
as the global density, the response and the correlation functions. On the other hand we have monitored some local 
quantities that allowed us to investigate the large-scale structures spontaneously emerging in these systems as a 
response to the imposed perturbation (driving). The comparison between global and local quantities allowed us to 
gain a deeper insight on how granular materials respond to perturbations and in this perspective of the importance 
of spatial structures. We have considered several different perturbation procedures defined in terms of the temporal 
functions describing the shaking amplitude (F in the experiments and x with F ~ l/log{l/y/x) in the models). In 
particular we have analyzed the case where one keeps x indefinitely constant and compared this case to the one where 
X varies as a function of time, x(t). Our main results can be summarized as follows. In the case of a procedure 
at constant x the system exhibits ageing described by a correlation function decaying as log{tw)/log{t + tw)- The 
response function exhibits a complex behaviour depending on x and ty^,. In general one always observes a positive 
response eventually followed, at times increasing with either x or by a negative response. All these properties can 
be explained looking at the heterogeneities arising in the density profiles. 

The scenario changes completely considering more complex shaking procedures where x — x{t). In this case (see 
Sect. IV for the details) the system can be found in an almost stationary state where ageing is not anymore present, 
i.e. the correlation function is time-translation invariant, and the response function is always positive. Also in this 
case the comparison of the results with the analysis of the density profiles allows us to gain a deeper understanding of 
the effect of the perturbation on the system. On this basis we are able to formulate some specific recipes for the best 
compaction procedure and to comment on some recent results concerning the validity of the Fluctuation-Dissipation 
Theorem and the possibility of a thermodynamic description for these non-thermal systems. 

Let us conclude with two points that open possibilities for future work: (i) recent experiments |]l^ have focused 
on the response, not to a slight change in the forcing, but rather to a large change, in the spirit of experiments of 



temperature cycling in spin-glasses 23 1. Work is in progress to reproduce the preliminary experimental results and 
draw comparisons with the spin-glasses phenomenology, (ii) we have defined the response function, as in as the 
response of the system to a change in the driving force, which is itself considered as an analog of the temperature. In 
order to try to understand if an extension of the FDT can be thought of, it would probably be more reasonable to 
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look at the response to a force acting randomly on the particles, and in a way uncorrelated to the overall driving, in 
a way similar to the one used in a lattice-gas model [|^ or in models of super-cooled liquids . 

Needless to say finally that it would be extremely important to have an experimental check of our predictions to 
be used as a starting point for further theoretical investigations. 
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